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Preface 

This report is a somewhat modified and expanded version of a report issued in June 
1985 by the School of Aeronautics and Astronautics at Purdue University entitled ’An 
Unsteady Lifting Surface Theory for Single Rotation Propellers’ [13]. The report 
provides background on a code, UPROP3S, that has been developed for the unsteady 
aerodynamic and aeroelastic analysis of advanced turboprops. Further information on the 
use of that code is given in a user’s manual [14]. The treatment of the blade boundary 
condition in Sec 3.1 has been altered from its original form, which was slightly in error. 
In addition, new sections on mistuning (Sec 1.4) and aeroelastic analysis (Appendix VI) 
have been added. The aerodynamic code based on this report has been incorporated into 
the aeroelastic analysis program ASTROP3 at NASA Lewis Research Center. However, 
the aeroelastic analysis given in App. VI is the one used in the code UPROP3S, and was 
not used in ASTROP3. Since the scheme has not been described elsewhere, a brief 
outline has been included herein. 
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Nomenclature 


C;j Aerodynamic influence coefficient, Eq.(28) 

Ci^ design lift coefficent of 16 series airfoils, Eq.(48) 
Op Thrust coefficient, Eq.(45) 

Cp Power coefficient, Eq.(46) 

D ± panel edge integrals, Eq.(30) 

Do part of D ± from Ko 
Dj part of D ± from Kj , Eq.(33) 

H Heaviside step function 

i %cr 

J advance ratio , ti/S 

K propeller kernel function, Eq.(II. 1 5) 

Ko quasi planar kernel function, Eq.(37) 

Ki residual kernel function , K - Ko 
— ) 

L normal to helical surface , Eq.(5) 

m interblade phase index, Eq.(20) 

M x axial Mach number 

Mo Helical Mach number of a blade point 
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n unit normal to blade surface 


Nr number of blades or periodic groups of blades 
P scaled pressure difference across surface, Eq.(24) 

Qij generalized aerodynamic force, Eq.(50) 
r h blade hub radius 

r t blade tip radius 

— > 

R blade surface rotation, Eq.(19) 

S tip to axial speed ratio, Qr t /U 
t time 

u perturbation velocity of fluid 

U axial flow velocity 

v n normal velocity of blade surface, Eq.(l 1 ) 

V side slip velocity 

W scaled upwash on blade surface, Eq.(24) 
x axial distance along axis of rotation, undisturbed air frame 
x axial coordinate in rotor frame, Eq.(2), Fig.l 
y cartesian nonrotating frame coordinate 
y cartesian rotating frame coordinate, Eq.(2), Fig.l 
z cartesian nonrotating frame coordinate 
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z cartesian rotating frame coordinate (along pitch change axis) 
Greek symbols 

a blade helical curvature , Eq.( 6 ) 

P angle between blade chord and plane of rotation 

P 3/4 blade setting angle at 3/4 tip radius 

A0 b angle between blades 

Ap front - back pressure difference 

8 blade surface displacement 

8 normal component of Eq.( 1 6) 

rj efficiency, Eq.(47) 

0 circumferential angle (nonrotating frame) 

0 circumferential angle (rotating frame) 

e control point position in fraction of panel chord 

% fraction of blade chord from leading edge 

po undisturbed air density 

o(r) blade surface helical coordinate, Eq.(4) 

co vibration, or excitation frequency 

co =co/Q 

Q. rotational velocity (positive in - 0 direction) 
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Section 1 
Theory 

(1.1) Introduction 

When a propeller undergoes blade vibration or operates in a nonuniform stream, 
unsteady aerodynamic forces are exerted on the blades. These fluctuating forces may 
influence the fatigue life, the aeroelastic stability and the noise output of the propeller. It 
is essential, therefore, to be able to predict aerodynamic responses for propellers in time 
dependent flows. This report describes a numerical method for making such predictions 
for single rotation propellers. 

The analytical model employed herein uses linear compressible small disturbance 
theory. In practice this means that the blades must be thin and at small local angles of 
attack (below stall). Further, any incident flow distortions must involve velocities which 
are small compared to the helical blade speed. Finally, linear theory does not account for 
embedded shocks which can occur on transonic tips or in a blade passage (at high 
solidity, where blockage may be significant). 

In addition to linearity, we assume that the load at any point on a blade fluctuates 
harmonically in time with some prescribed frequency 0). Hence the entire disturbance 
field will fluctuate harmonically at the same frequency (in a frame rotating with the 
propeller). More complex periodic disturbances can, of course, be examined by breaking 
the disturbance into its Fourier components and finding the response to each individually. 

A blade may be thought of as a mean camber surface overlaid with a thickness 
distribution. In the linear approximation the thickness distribution does not effect the 
loads and will, therefore, be ignored. The mean camber surface, which determines the 
loads, may be either rigid or may vibrate harmonically about some average position. In 
either case the object of the calculation is to find that load distribution (pressure jump 
across the surface) for which there is no flow through the surface. 

The task of finding such a load distribution is simplified by transferring the 
boundary conditions from the actual camber surface to a neighboring helical surface (see 
Sec. 1 .3). This is done by choosing some curve, called the "generator," on the surface. If 
the blade vibrates, the generator is chosen to lie on the average camber surface. As the 
blade advances and rotates the generator sweeps out a fixed helical surface. If the actual 
blade lay exactly on the helical surface - and there were no incident flow nonuniformities 
- then the blade would produce no disturbance as it slices through the air. Thus it is the 
deviation of the camber surface from the helical surface swept out by the generator which 
produces a load. 

In practice, the load distribution is placed on the helical surface, rather than on the 
camber surface. The advantage of this transfer is that the trajectories of all points at a 
fixed radius on the helical surface are the same. Since, as will be seen later, the 
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calculation of the velocity induced by one load point on the blade requires an integration 
along the trajectory of that point, the present scheme allows one integration to serve all 
blade points at a fixed radius. If the load distribution were placed on the camber surface 
then a separate trajectory integral would need to be done for each blade point. 

In this model, then, a propeller blade is thought of as a distributed load sliding 
along the helical surface which induces a given distribution of normal velocities on itself 
(Sec. 1.3). 

Linear aerodynamic theory provides an explicit expression for the normal velocity 
at any point in terms of a surface integral of the load distribution. Since the velocity is 
given, this relation amounts to a linear integral equation for the load (Sec. 1.5). To solve 
this integral equation the blade is broken into a finite number of elements, on each of 
which the load is considered to be constant. To each element a control point is assigned 
at which the normal velocity is specified. The loads on all the elements are then 
determined simultaneously by requiring that their resultant induced normal velocity have 
the specified value at each control point. 

The principle difficulty in implementing this scheme is that the velocity field 
induced by a constant load blade element cannot be computed exactly. In fact, not even 
the velocity induced by a point load can be computed exactly (meaning, of course, that 
evaluating either quantity to machine accuracy would require excessive computation 
times). Hence the influence coefficients (the normal velocity at a control point induced 
by a unit load on an element), must be approximated numerically with enough accuracy 
to make the results meaningful yet with enough simplicity to make the scheme practical. 
The way in which the influence coefficients are calculated is the key to making the 
method work. The scheme used for the present calculations is described in detail in 
Chapter 1.6. 

Several three dimensional linear lifting surface theories for propellers, similar in 
various respects to the present method, have appeared in the past. The work of 
Hammond et al. fl] is closest to the present scheme, differing only in the discretization 
procedure. Hanson [2] uses the same theoretical framework but expresses the kernel 
function as an infinite series of Bessel functions. He gives the equations, but no results, 
for a vibrating blade. Farassat [3] and Long [4], employ linear compressible theory, as 
does the present work, but do not use the mean surface approximation. Like Hanson, 
they show no results for a vibrating blade. Sullivan and co-workers [5-6] have developed 
a vortex lattice method for steady performance calculations (recently extended to 
counter-rotation, using a quasi-steady approximation [7]. Finally we note that a number 
of 3-D finite difference algorithms have been developed for the steady problem [8], 
though none has been extended to unsteady flows. 
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(1.2) Coordinates and Geometry 


Coordinate Systems 
(see Fig. (1)) 

(a) Inertial Frame 

(x,y,z) is a right handed cartesian system fixed in the undisturbed fluid. In this 
frame the propeller center of rotation advances with speed U along the (-x) axis and the 
blades rotate, clockwise, at speed Q when viewed from the -x axis (so the rotational 
velocity vector points along the +x axis). 

(x, r, 0) are cylindrical coordinates, with 0 measured from the z-axis; so that: 

y = rsin0, z = rcos0 (1) 


(b) Rotor Frame 

(x, y, z) is a right handed cartesian system fixed in the rotor. (The z axis can be 
taken to coincide with the pitch change axis of the reference blade). 

(x, r, 0) are cylindrical coordinates fixed in the rotor. 

From the sign convention adopted we have, 

x = x + Ut , 0 = 0 + Qt 

y = r sin 0 , z = r cos 0 (2) 


Helical Surface 
(see Fig. (2)) 

Any fixed point on a blade follows a helical trajectory x-U/Q0 = x- U/ Q0. 
Let x g (r), 0 g (r) define some space curve lying on the blade chord surface ( e.g. 
The trailing edge or mid-chord line). This curve generates a helical surface, 

x = U / Q (0 + a) (3) 

where 

a(r) = Q / U x g - 0 g . (4) 

Note that this helical surface is time-independent, and fixed by the advance ratio 
and the choice of generator curve (which determines o(r)). 

The normal to the helical surface is, 
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L = V(x - U / Q (0 + a)) 


=T X - U / (Qr) 0 9 + od r ) 


( 5 ) 


where (i x , ie, i r ) are unit vectors in the (x, 0, r) directions respectively, and 


a = r 


do 

dr 


( 6 ) 


Note that the normal L is a function of position only. 

We shall parameterize the helical surface by (r, 0). An area element on the surface 
is, then, 


dA = L rdr d0 


(7) 


and an element of arc length along the surface at fixed radius is, 

ds = d0 [r 2 + (U/Q) 2 ]- 5 


( 8 ) 


In the aerodynamic model, the load distribution is placed on the helical surface, 
within the region 


r h < r < r t 

®LE( r ) ^ 0 ^ QteOO • 


(9) 


These edge coordinates of the lifting surface are obtained by projecting the blade 
chordline onto the helical surface so that the (helical) arc lengths are equal: 


0TE ~ 0LE = 


(Ax) 2 +(rA0) 2 

(U/fi) 2 +r 2 


(10) 
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(1.3) Surface Boundary Condition 

The normal velocity of the fluid at the blade surface must equal the normal 
velocity of the blade surface at every point. Let u be the fluid velocity in the inertial 
frame ( i.e. relative to the undisturbed fluid), 7? the unit normal to the blade surface, and 
v n the normal velocity of the blade surface (again in the inertial frame). Then at every 
point on the blade we must have, 

Tvu = v n . (71) 

In the present (linearized) model, wherein the load is transferred to the helical 
surface, the boundary condition is similarly transferred. Thus at every point on the 
helical surface we impose the constraint, 

LTi= |L|v n , 02) 

which says that the normal velocity induced on the helical surface by the load 
distribution must equal v n (which is known once the blades shape and motion are 
specified). 

In the present study the normal velocity is assumed to be either steady or simple 
harmonic at any fixed blade point. Expressions for three cases will be given here: steady 
operation, propeller in yaw, and blade vibration. 

Steady Operation 

For steady state operation we have, 

v n = -Un x - Q(z n y - yn z ) , (13) 

(where n x , n y , n 7 are the components of n* in the blade coordinate system). 

Blade in Yaw 

If the propeller blades are rigid but operate in a crosswind V (in the y direction), 
the corresponding amplitude of v n (on the reference blade) is: 

v n = -V(n y + in z ). (14) 

The fluctuating frequency is co = Q and the interblade phase index is m = -1 (see 
"periodicity".) This case is equivalent to a yaw angle of tan -1 V / U = V / U. Of course, 
the yaw angle must be small in order for linear theory to be applicable, in which case the 
resultant blade loads will be proportional to the yaw angle. 

Blade Vibration 

If the blades vibrate with frequency to, the amplitude of the resulting normal 

— ) 

velocity fluctuation may be computed from the displacement vector 8 and rotation vector 
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R assigned to each blade point!?: 

v n = ico5 + v R (15) 

where 

8 = nf> (16) 

is the normal displacement, and vr = v R[ + vr 2 , with 

V R, — U ( n y R z n z Ry ) (17) 

+ Q[(z n z + y n y )R x - n x (y R y + z R z )] 

being the contribution from rotational deformation, and 

— £1 ( 8 y n z 8 Z n y ) (18) 

being the contribution from the increase in speed due to a radial displacement (which was 
missing in the original report and code.) 

These expressions can be used directly if the structural displacement and rotation 
vectors are available (e.g. from a finite element analysis). If the rotation vectors are not 
available (or are unreliable), they can be computed from their fundamental definition in 
terms of derivatives of 1), 

R x = -n y n't 7 + n z ivf> y ( 1 9) 

= n x n*^ z 
R z = -n x m8 y 

The subscripts on 8 denote differentiation along the indicated directions, assuming the 
blade is parameterized by y and z.(This method of computing rotation vectors was not 
used in the original work and is superior to the scheme described earlier.) 

In practice, the blade is usually defined by a discrete set of nodes covering the 
camber surface, with !?, 8, and .perhaps, R defined at each node. To determine the 
normal velocity at any prescribed point on the surface, the!? data is fit with a polynomial 
( quadratic in y and bilinear in z and y) .least squares interpolated to the six nearest 
neighbor nodes. The five coefficients of this surface fit can be used to define the surface 
normals as well as the interpolations of the displacement vector and its derivatives. 
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(1.4) Periodicity 

We assume that there are Ng identical blades (j=l, Ng) spaced equally in angle 
increments A0g = 2rc/Ng. Then the load on any one blade will be the same as the load on 
every other blade except for possible phase shifts. To allow for such phase shifts we 
assume that the normal velocity on the jth blade at a given point is related to the normal 
velocity on the reference blade (j = 1) at the corresponding point by, 

(v„)j = (v„x e™ 0 -'’ 48 - (20) 

where m may assume any integer value, m = 0, 1, ..., Ng— 1. 

Let Ap be the pressure difference (upwind-downwind) across a blade. Then, by 
(17), the loads are also related by 

(Ap)j = (Ap)] e im(H)A0B , (21) 

so that only the load on the reference blade need be found. 

Note that with Ng blades there will be Ng helical surfaces with helix numbers 
(Eq.(4)) given by, 

Oj = o, -(j-1) A0 b (22) 

Finally we shall assume that the normal velocity and load at any fixed point on a 
blade vary harmonically in time, with frequency 0), like e ,<Bl . For most interference 
problems co will be a harmonic of the rotor frequence fi. For vibration problems, co will 
be close to a blade natural frequency. 

Mis tuning 

The foregoing discussion is easily extendible to cases in which not all the blades 
are identical, but the ’mistuning’ is distributed periodically( eg. alternate mistuning.) 
This is done simply by redefining ’blade’ to mean ’a group of blades’. The groups are 
assumed to be indistinguishable, but each group may consists of an arbitrary collection of 
disjoint surfaces. (The only special coding change that was needed to allow for different 
blades within a group was to ensure that the tip of one blade was not connected by panels 
to the root of the next blade.) If all the blades are different (i.e. there is no periodicity) 
then there is ,by definition, only one group. A ’tuned’ rotor consists of Ng groups with 
one blade per group. Alternate mistuning would have Ng groups with two blades per 
group. 

In addition, mistuning studies require a generalization of ’blade mode’. The 
simplest definition of a ’group mode’ is to suppose that one group mode corresponds to 
one blade in the group vibrating in one of its natural modes while all other blades in the 
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group remain fixed. Thus the total number of group modes will be the number of modes 
retained per blade summed over the group. The numbering sequence is arbitrary. 
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(1.5) Lifting Surface Integral Equation 

From the assumptions of aerodynamic linearity, simple harmonic time dependence, 
and blade to blade periodicity, the entire disturbance field surrounding the propeller can 
be determined directly from the distribution of load amplitude over the reference blade. 
In particular, then, the load distribution determines the distribution of normal velocity 
over the reference blade. 

This relationship between load and normal velocity can be expressed as an integral 
equation, 


1 ^TE 

W(r, 0) = J j P(r c , 0 o )4r K((M) 0 , r, r o )r o d0 o dr o 

r h 9| .f. 


(23) 


where 

7 * — * 

W = 4rc(“^y~)e iw0 

p = S 2^E_e i “ §<> . (24) 

PoU 2 

Here W is proportional to the normal velocity, P to the pressure jump across the 
blade, and K is a (known) kernel function. In most applications the normal velocity is 
specified and the integral equation must then be solved for the load distribution. The 
method used to do this will be described in subsequent sections. 

Expressions for the kernel function are derived in Appendices I and II. Here we 
will only discuss its general properties. 

Parametrically, K depends on the speed ratio, S; the frequency ratio, co = co/Q; the 
axial Mach number, M x ; the interblade phase index, m; the number of blades, N B ; and 
the helix numbers, a(r). If these quantities are specified then K is uniquely determined. 
The kernel does not depend on the blade planform shape or the normal velocity 
distribution. 

Physically, K can be interpreted as follows. Suppose that the load is concentrated 
on a particular radius, with zero amplitude upstream and constant amplitude downstream: 

P(r, 0) = H(0)5(r-r o )/r o . (25) 

If there are N B blades then there are N B such loaded lines, all extending from the same 
axial position to downstream infinity, with equal strengths but phase shifted in accord 
with the index, m. The normal velocity induced by this set of loaded lines on the 
reference blade is, from Eq’s. (23) and (25), 
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W(r, 0) = K(0, r, r 0 ) . 


(26) 


The kernel depends only on the angle difference 0-0 o , because of rotational 
symmetry. It depends on the radii r, r 0 separately, because the relative velocity is a 
function of radius. The evaluation of K as a function of 0 - 0 O for any specified pair of 
radii r, r G , requires numerical integration along the Nb helical lines at radius r„ (the 
integration is a consequence of Newton’s law - velocity is the integral of the force 
producing it). As a result the evaluation of K tends to be the dominant part of the load 
calculation for moderate panel density. When the number of panels is very large the 
0(n 3 ) cost of Gaussian elimination will dominate. 

From the physical interpretation of K as the velocity induced by a set of line loads 
it is apparent that K will contain singularities. In particular, K will always be infinite on 
r = r 0 , 0 > 0 O , i.e. on the loaded line. If the helical Mach number at radius r 0 is less than 
one then K will be bounded everywhere else. However, if the helical Mach number is 
greater than one, then a spiral Mach cone will extend back from the tip of each of the 
loaded lines. Thus K will be infinite wherever one of these Mach cones intersects the 
reference helical surface. The singularities of K must be handled very carefully in 
constructing a numerical solution of the lifting surface integral equation. 
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(1.6) Discretization of the Integral Equation 

In brief, the lifting surface integral equation, Eq(23), is solved approximately by 
splitting the blade into a finite number of elements within each of which P is assumed 
constant. The normal velocity W is then specified at one point per element, thereby 
reducing the integral equation to a set of simultaneous algebraic equations for the loads 
on each element. 

More specifically, the blade is split into NRP radial strips of arbitrary width. Each 
strip is then divided into NXP chordwise pieces by a sequence of constant partial chord 
lines (the chordwise spacing is arbitrary). The blade is thereby decomposed into 
NP = NRP-NXP quadrilateral panels. 

A control point is placed on each panel at its mean radius, a distance eA0 
downstream from the midspan panel leading edge, ( where A0 is the panel chord at 
midspan , and distance is measured in 0.) The arrangement is illustrated in Fig. (3). The 
selection of e is discussed in Sec. (2.2). 

The algebraic system resulting from this discretization is 

NP 

Wi = I CijPj (27) 

j=i 

where 


W; = W at ith control point, (0i, r;) 


Pj = P on jth panel, 

9K(0 ; -0o.ri.ro) 

Qi = -jj d0 o r o dr o 

00o 


(28) 


Note that the chordwise integration in the influence coefficient C;j can be performed 
directly (this is why P, rather than Ap, was treated as being constant over the panel). 
Hence 


Qj = D_ - D + (29) 

where D ± are the radial integrals along the panel leading (-) and trailing (+) edges 

D ± = J K(0i - 0 O± , ri , r 0 )r 0 dr 0 (30) 

(0o±(r o ) are the panel leading and trailing edge coordinates). Note that for a fixed control 
point, D + for one panel is the same as D_ for the next panel in the row. Thus D must be 
evaluated NXP+1 times per radial strip. 
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(1.6.1) Numerical Evaluation of Influence Coefficients 

Computing an influence coefficient requires a radial integration of the kernel 
function for D. However, the kernel function evaluation itself must be done numerically 
so that only a few values of K will be available near any given panel. This is not a large 
problem, if the panels are small, unless the panel encloses, or lies near a singularity of the 
kernel. This does happen, however, when the control point radius lies within, or near, the 
panel row. 

To circumvent at least some of the difficulties associated with the kernel function 
singularity, the kernel is split into two parts 

K = K 0 + Kj (31) 

with a corresponding decomposition of D: 

D = D 0 +D,. (32) 

The Kq term is an analytical approximation which is valid when r-r„ is small, and 
therefore contains the dominant part of the singular structure of K (to be precise, K and 
K<, behave like (r-r 0 ) -2 near r = r 0 when 0>9 O ). The residual kernel, Kj, is simply 
defined as K-Kq. For the cases run, Kj has been found to be also singular, but only as 
(r— r D ) -1 . This singularity will be accounted for in the evaluation of Di . 

The specific form of K<, used here will be given in the next section. The significant 
point is that it is simple enough that D 0 can be evaluated exactly and efficiently. 

The part of the influence coefficient arising from the residual kernel K] , i.e., Dj , is 
evaluated by simple quadrature rules. The specific rule used depends on whether the 
control point lies within the panel row or not. (Cases (a) and (b) below) 

Let the upper and lower radii of the panel be R] and R 2 respectively and suppose 
the particular panel edge in question runs from 0 O ] to 0 o2 (see Fig. 3). Then the value of 
Dj required is: 

^2 

Dt = J K, (0 - 0 O , r, r 0 )r 0 dr 0 . (33) 

R. 

(a) r*(R 1+ R 2)/2 = R m 

In this case the control point is outside the panel row. The kernel K] is evaluated 
numerically along the midspan and Dj is evaluated from 


D, = R m (R 2 - R, )K, (0 - 0 om , r, R m ) Q 


(34) 



where R m and 0 om are the mid-points as shown in Fig. 3. If 0< 0„ m , ie. if the control 
point is upstream from the panel edge, then the factor Q is set to 1. If 0 > 0 om , though, 
the residual kernel K! has a pole at 7 0 =7 which may be close by. In this case we take 


Q = 


R m ~r 

R 2 -Ri 



r, r„-r 


Rm-r 

R2-R1 


In 


R2 — r 
Ri-r 


(35) 


If the control point radius is more than one or two panel spans removed, this correction 
factor is nearly one. It differs from one significantly only in the neighborhood; e.g. if 
7= R m + R 2 - R] then Q = 1.1. 


(b) 7= R m 

In this case the control point is in the panel row and the kernel is singular along the 
midspan. The kernel is then evaluated numerically along the edges R] and R 2 and D] is 
evaluated from the quadrature rule: 

D, =(R 2 -R,) y R m [K!(0-0 o2 ,7, R 2 ) + K,(0-0 o 1 ,7, R,)l (36) 

Note that this rule is exact if K! = a/(r 0 -r) (i.e. it gives Dj =0) and if Kj is linear 

in7 c . 

In summary, if the control point lies outside the panel row we evaluate the residual 
kernel along the midspan of the row. If the control point lies within the row then we 
evaluate Kj along the outer edges of the row. For each control point and panel row we 
will need numerical values for K! at a discrete set of angles 0-0 o , numbering roughly 
(NXP) 2 . In practice a single table of Kj values covering the needed range of 0-0 o is 
generated for each panel and control point row and values for specific panels and control 
points are found by interpolation in the table. 
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(1.6.2) The Quasi-Planar Kernel, K<, 

As observed earlier the kernel K is the normal velocity on the reference helical 
surface induced by a system of loaded lines emanating from each blade. If the observer 
is very close to the loaded line on the reference surface, the influence of the other blades 
will be negligible. Moreover the effect of curvature will be small if the distance from the 
line to the observer is small compared to the radius of curvature of the helix. Hence the 
local effect should be essentially the same as that of a straight loaded line on the tangent 
plane. Thus we shall use a modified form of the planar wing kernel function. 

Let M 0 = [U 2 +(Qr 0 ) 2 ] v Va o be the helical Mach number at radius r. The form of 
Kq depends on whether M 0 is greater or less than one: 


Ko(A0,r,r o ) = 



{X + [X 2 +BY 2 ]} sup half if M 0 < 1 


= 2CH(X) [X 2 +BY 2 ] l/ * if Mo > 1 

■=■ v2 


(37) 


where 


B = 1 - M 2 


M, 


Mx \2 1 


C =— 

m o r 


M„S 


(38) 


and H is the unit step function, and a is the helical surface curvature at radius r defined 
by Eq. (6). 

If we interpret X as being distance downstream and Y perpendicular to a uniform 
flow at Mach number M 0 , then these are the planar wing kernel functions for steady 
subsonic and supersonic conditions respectively (the supersonic kernel is taken to be zero 
outside the Mach cone). 

In terms of propeller coordinates we take 

X = A0 + (M X /MJ 2 a(Ar/r) 

Y 2 = (M,/M 0 ) 2 [1 + (M x /M 0 ) 2 a 2 l(Ar) 2 (39) 

Ar = r-r 0 

It should be noted that the precise form of Kq is rather arbitrary. The only 
requirement is that K 0 should capture the dominant singularity of K and be integrable 
analytically. The present choice has these properties, though there may well be better 
choices available. 


14 



The corresponding values of D 0 can readily be evaluated by use of the following 
identities. Let A9 be linearly interpolated along the panel edge so that X = X 0 + X Y 
where X 0 and X' are interpolation coefficients determined from the panel edge end 
points. Let R(Y) = [X 2 +(1-M 2 )Y 2 ] 1/4 . Note that Y is linear in r 0 so that integrating in 
r 0 is equivalent to integrating in Y. Then, if d = (X ) +B: 

*111 = {_2^._x'ln^i?-+dln(X , X+dR+BY)} (40) 

Y 2 dY v Y Y 2 


gives D 0 in the subsonic case; 


JL = -iL {- — -X'ln *±!L + d In (X'X+BY+dR)} 
Y 2 dY 1 Y Y 


(41) 


gives D 0 in the supersonic case (subsonic edge, d > 0); and 

R d , R X+R . (J1 . , X'X+BY N1 

— — = { X In — — + d sin ( nr)) 

O 1 a/ * 1 — ' M 


dY 


[-BX 2 ]* 


(42) 


2 

gives D 0 in the supersonic case (supersonic edge, d < 0). 
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Section 2 
RESULTS 

(2.1) Introduction 

The theory outlined in Section 1 of this report has been implemented in a program 
UPROP3S (unsteady-propeller- 3D-single rotation). Representative results based on that 
program will be described here, with a view toward demonstrating and verifying the 
program’s capabilities. 

In outline the cases which will be discussed are: 

1 . Convergence Study 

The method is shown to converge when the number of panels is 
increased. A control point location for best convergence rate is 
identified. 

2. Planar Wing 

The method is shown to reproduce known results for isolated planar 
wings in steady and unsteady flow. 

3. Performance Characteristics 

The theory predicts steady state propeller performance 
characteristics which are in agreement with existing methods. 

4. Ten Bladed Fan Study 

The analysis of a 10 bladed fan shows the same qualitative 
influence of interblade phase angle on the vibratory response 
loads as 2-D cascade theory. 

5. SR3 propeller 

Generalized forces for structural vibration of an advanced 
turboprop design are presented. 
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(2.2) Convergence and Control Point Position 

The lifting surface integral equation, Eq(23), admits an infinite number of 
solutions. However, only one is physically acceptable - the one which satisfies the Kutta 
condition, Ap = 0 at the trailing edge. Since this condition is not clearly imposed in the 
numerical scheme, the ability of the method to capture the Kutta condition needs to be 
demonstrated. To be precise, we must ask the question: does the discrete solution 
converge to the desired solution of the continuous problem as the number of panels is 
increased? 

Numerical experimentation has indicated that the answer to this question is yes - 


provided that e > y (i.e. as long as the control point is downstream from the panel 
midpoint). Thus the numerical Kutta condition is simply £ > — . (This problem is further 


discussed in Appendix V.) 

We illustrate the convergence properties of the method with an oscillating 
propeller problem: 


Purdue model blade (Table 1) 


N b = 1, S = 1.7, M x = .1, (0= 1 

v n = Qr t (l+ico0) (43) 

A sequence of calculations were performed in which the number of radial panel rows was 
fixed at 9 and the number of chordwise panels (NXP) and the control point position (e) 
were varied. 

Fig. (4) shows the fluctuating sectional thrust coefficient (real and imaginary parts) 
at 80% tip radius. For each value of e, results were calculated at NXP = 4, 8 and 1 2. The 
data was then fit with a quadratic in 1/NXP, as shown. The zero intercept is a prediction 
of the limit value. 

The important conclusions to be drawn from this plot are: 

• The limit values are insensitive to control point position. 

• The choice £ = .85 gives nearly the same result for all NXP. 

• The algorithm is first order (error proportional to 1/NXP.) 

Thus the numerical scheme is convergent and the optimal convergence rate is obtained 
with £ = .85. Similar results for the 2D airfoil problem are given in Appendix V. 

The in-phase and out-of-phase parts of the chordwise pressure distribution at 7= .8 
are shown in Fig. (5) for the £ = .85 case. It is clear that the Kutta condition is captured. 


17 


Moreover the pressure distribution (as well as the thrust) is insensitive to NXP at this e, 
though of course the larger NXP the better defined the distribution. 

All other results reported here will be for e = .85. 
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(2.3) Planar Wing 

The propeller analysis can be used to approximate an isolated planar wing either 
by taking the speed ratio S (Q r t /U) to be very small or by taking the blade span to be 
small compared to the radius of rotation. The results of such a simulation are shown in 
Figs. (6, 7) with comparisons to Albano and Rodden’s doublet lattice method. [91 The 
case chosen is an aspect ratio 2 rectangular wing at low Mach number, pitching 
harmonically about midchord. (Simulated so that the variation in flow speed along the 
span is less than 

Figure (6) shows the real and imaginary parts of the chordwise load distribution at 
midspan for a reduced frequency coc/U = 1. Fig. (7) shows the total lift coefficient per 
unit amplitude as a function of reduced frequency. 

It is apparent from these comparisons that the present analysis does produce the 
same results as classical linear unsteady wing theory in this limit. 
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(2.4) Propeller Performance 

In this section we will discuss the steady state performance characteristics of three 
propellers, comparing the predictions of the present lifting surface theory, and Sullivan’s 
vortex lattice method [5-8] which is based on linear incompressible (steady) 
aerodynamics and so should agree with the present theory at low Mach numbers. Any 
difference between the two theoretical predictions (at low Mach numbers) then, is a 
consequence of the different discretization methods. In this regard it should be noted that 
the two theories differ significantly in the way efficiency is calculated: the present 
scheme uses the leading edge thrust concept while Sullivan [5] uses the Kutta-Joukowski 
law. 

Results will be presented in terms of the conventional performance parameters: 


AdvanceRatio : J = it/S 

(44) 

thrust coefficient : Cj = thrust • 7t 2 /(4Q 2 rf ) 

(45) 

power coefficient : Cp = power • rc 3 /(4Q 3 r?) 

(46) 

efficiency : T] = JGr/C p 

(47) 


The three configurations studied (all straight bladed propellers), are: 

(a) Purdue model (Table 1) , Nr = 1 

(b) NACA 109622 (Table 2), N B = 3 

(c) SR2 (Table 3), N B = 8. 

The chord, twist and camber distributions are given in the indicated Tables. (AP is 
the change in blade setting angle from the 3/4 radius station.) 

The normal velocity used in the lifting surface calculation is, 

v n =U r [p- tan- 1 (■£-) + -^ In (-^-)1 . (48) 

i2r 4tt 1-q 

where £ is the distance along the chord line measured from the leading edge, in units of 
chord. The term in this expression is based on the camber line of the 16-series 
airfoils used on the outboard sections of the SR2 blade (the present calculations use the 
same airfoil at all radii). The value assigned to is the ’design’ lift coefficient of the 
airfoil at zero angle of attack (the 16 series airfoils being designed to give a flat load 
distribution according to linear 2D airfoil theory.) The first two blades are uncambered 
(CLp = 0) and so derive their thrust from twist only. The SR2 blade is cambered and so 
derives part of its thrust from the camber. 
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Figure ( 8 ) shows the power and efficiency as a function of advance ratio for the 
Purdue blade operating at low Mach numbers and P 3/4 = 45.7°. The two methods predict 
essentially identical efficiencies. The present method gives a power coefficient slightly 
(around 10%) lower than the vortex lattice method. (The thrust, which is not shown, 

must also be slightly lower to get the same efficiency). 

Figure (9) shows a similar comparison between the present method and vortex 
lattice theory, for the three bladed NACA 109622 (again at low Mach number). The two 
results are clearly very close to each other, with the lifting surface theory now slightly 
lower in terms of efficiency as well as power and thrust. 

The general conclusion to be drawn for Figs. ( 8 ) and (9) is that for straight bladed 
propellers operating at low Mach numbers the present theory predicts essentially the 
same performance characteristics as the vortex lattice method. 

In Fig. (10) we compare the power coefficient predicted by the present scheme to 
measured values for the SR2 propeller. Three configurations are shown: 
(M x , p 3/4 ) = (.3, 40°), (.3, 52°) and (.7, 58°). In each case the relative tip Mach number 
at the lowest advance ratio is near 1. Thus compressibility effects are to be expected. 
For each blade setting two theoretical predictions are shown - one at low Mach number 
and one at the experimental value. As would be expected, the incompressible result is 
low. The predictions using the appropriate Mach number are in very good agreement 
with experiments. 

This indicates that the linear compressibility effects included in the present theory 
are a significant factor. 
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(2.5) Torsional Response of a JO Bladed Fan 

In this example we consider a fan with 10 blades, each executing a torsional 
oscillation. The objective of the study is to validate the present theory with regard to the 
dependence of vibratory response on interblade phase angle. 

The blade planform geometry is given in Table 1. The blades are unswept, 
constant chord, and helically twisted. The operating conditions are defined by M x =0.1, S 
= 1.597, co = 1.0. 

The normal velocity distribution on the reference blade is, 

v n = U re i .a(l+ia>0) (49) 

where U re i is the helical velocity and a is the local angle of attack amplitude, which is 
assumed to vary linearly from hub to tip a =7-7 h • 

The aerodynamic response, shown in Figs. (11,12), is measured by the sectional 
lift coefficient l/jtpb 3 0) 2 , where b is the blade semi-chord (.23 r t ,) and 1 is the lift per unit 
span. Results from Smith’s two dimensional cascade theory [10] are included for 
comparison. 

The two figures are for two radii: one inboard, T = .6 and one near the tip, 7 = .9. 
In both cases the 2 and 3 dimensional theories show the same qualitative dependence of 
load on interblade phase angle. Not surprisingly, the agreement is quantitatively better 
away from the tip. The 2-D analysis drastically over predicts the inphase loads at 7 = .9 
for most interblade phase angles. This is to be expected for an unducted fan and is one of 
the reasons why a 3-D analysis is preferable. 
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(2.6) Generalized Forces for the SR3 

In an aeroelastic calculation the aerodynamic loads induced by blade vibration 
appear as generalized forces, 

Qjk = j j A Pk 5j dA (50) 

where 8j is the normal displacement amplitude for the j* mode and Ap^ is the amplitude 
of the pressure difference across the blade resulting from motion in the k 111 mode. The 
vibration is assumed to be simple harmonic with frequency to. 

Generalized force results will be presented here for the SR3CX2, a swept, flexible, 
advanced turboprop blade designed for flutter testing. In-vacuum vibration mode shapes 
(including the effects of centrifugal stiffening) were obtained from MSC NASTRAN for 
6100 RPM rotational speed with p 3 / 4 = 60.7°. These results were provided to the author 
by NASA Lewis Research Center. 

The normal displacements in mode one are shown in Fig. (13). The generalized 
force Qu is shown in Fig. (14) as a function of frequency. 
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Appendix I 

Velocity Induced by a Moving Point Force 


Consider a uniform unbounded stationary fluid with density p G and sound speed a 0 . 
We imagine that this fluid is disturbed by the action of a point force moving along some 
prescribed path. The resulting pressure disturbance, p', and velocity, u, are presumed to 
obey the linear acoustic equations of continuity and momentum conservation: 

-^l + p o aJV-5=0 (U) 



+ Vp' = F(t) 5(x-x 0 (t)) 


( 1 . 2 ) 


where -F is the force exerted on the fluid and x 0 is its point of application. 

The solution of these equations for p' has been given by Lowson [11]. In this 
application, though, we need the velocity. To get it we introduce a "displacement 
potential," \| f (such that V \| t is the particle displacement): 

u = V\|/ +u s (1.3) 


P = -Po 



(1.4) 


where is an impulsive velocity concentrated on the trajectory. Substituting these 
definitions into the momentum equation, Eq(1.2) yields an expression for u^. 

= — F(t) 8(x-x 0 (t)) (1.5) 

3t po 


while, from the continuity relation (1.1) we get: 


d_ 

at 


[V 2 V- 


1 3 2 V,_ 


a o 


at 2 


1 = V.U 5 


(1.6) 


Integrating (1.5 and 1.6) over time yields: 

oo 

u 5 = — J FdOfiCx-XodOJHO-t^t, 

Po -oo 

2 00 

v 2 y- -V = — v - 1 F(ti )6(x-x 0 (ti ))(t-t] )H(t-t, )d tl 

ao at 2 Po _oo 


(1.7) 

( 1 . 8 ) 
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where H is a Heaviside step function. 

Note that U 5 vanishes everywhere except on the trajectory of the point force. 

The solution of (1.8) for the displacement potential can be obtained from the point 
source solution to the wave equation: 



9 f(t-ti-Ri/ao) 
3t 2 } 4 tcR, 


= -5(X-X 0 (t! ))f(t t ! ) 


(1.9) 


where R, = | x-x 0 (t, ) | and f(t) is an arbitrary function. From this it is apparent that y is 
given by: 


1 

\u = 

47tpo 


oo 

V. J F(t,) 


(t-ti -Ri /a<, )H(t-tj -Ri /a 0 ) 


dti 


( 1 - 10 ) 


After performing the divergence operation we get, 

— > — > 

00 F 

\|T= -T- — I — (t-t 1 )H(t-t 1 -R 1 /a 0 )dt] (1.11) 

4tup 0 jL R] 

The step function in the integrand means that we need integrate only over that part of the 
trajectory which is acoustically accessible to the field point. We may remove the step 
function by introducing the retarded time x, which is the time at which a sound wave 
must be emitted from the source in order to arrive at the field point x* at time t. For fixed 
field point position and time, the retarded time must be determined by finding the roots of 
the transcendental equation, 

t = T + R(T)/ao. (1.12) 

Differentiating this relation with respect to x yields 

A = Jl = i-M R (M3) 

dx a D R 


where Tt 0 =3t 0 (x) is the velocity and Mr is the Mach number in the direction of the 
observer. Clearly if the force moves subsonically then there can be only one x for each t 
(since Mr is then everywhere less than one). If the source Mach number is greater than 
one, though, there may be several retarded times (or none) which satisfy Eq(1.12). 

In general, let Xj, X 2 , X 3 ... etc. be the ordered roots of (1.12). Let us further 
assume that Mr <1 at Xi. Then we must have Mr > 1 at X 2 , Mr < 1 at X 3 etc., since the 

extrema (Mr = 1 or — = 0) must occur in between roots. It follows that for any 

dx 

function f(t): 
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oo Tj ^5 

J fOi )H(t-t|-Ri/ao)dti = J f(ti)dtj + Jf(tj)dti + Jf(tj)dti + 


= / f(tj )dt t - J f(tj )dt! + Jf(t,)dt! 

— oo — oo — oo 

1_M R f r/ x, 

= ?Ti^rl f(,|)d " 

where Mr is evaluated at the root x and the summation is over all possible roots. 
The result (1.1 1) for the displacement potential can thus be expressed as: 

1 1-M R ] 

To get the particle displacement we need the result (from (1.12)) that 

— ) 

*7 N 

aod-M R ) 

N = R/R 

where N is a unit vector directed from the source at x 0 (x) to the observer. 

It follows directly that the velocity is 


-> 3 f 1 

u = — [ 


1-M r 

I TT-TT^ G(X)] 


at 47tp G 7 |i-Mr| 


0.14) 


(1.15) 


(1.16) 


(1.17) 


where 


G(x) = - 


N(F-N) 
a*R(l-M R ) 


+ J -rf [Ft^NiCFrNOldt, 


(1.18) 


— 7 v 

Note that G (and therefore u) is infinite on the trajectory (since Rj = 0 on the path 
of integration), and, in the supersonic case, on the Mach cone Mr = 1 . 

Equations (1.17-18) give the fluid velocity induced by a varying point force moving 
along an arbitrary path. 
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Appendix II 

Kernel Function for Helical Lifting Surfaces 


The results of Appendix I are for an arbitrary time dependent force moving along 
an arbitrary path. We now specialize to the case of helical motion and a harmonically 
fluctuating force. The path of the force, is, then, 

x„(t) = x c -Ut , e o (t) = 0 o -Qt, r G = const. . (II.l) 

The force, F is perpendicular to the helical surface at x^, so if Ap(r 0 , 0 O ) is the amplitude 
of the pressure difference across the surface, then 

F(t) = fLo e iox 


where, 


f = Ap(r 0 , 0 O ) r o dr o d0 o 

L 0 =L(x c (t)) 


(H.2) 


(L is the surface normal defined in Eq. (5).) 

Now let L = L(x) be the normal to the helical surface at the field point x. From 
Eq.’s (1.17,18) and (II. 2), the induced normal velocity is 




1-M r 
I 1-Mr | 


(II.3) 


where 


G = - 


(LN)(L 0 N) + j H, 

a oR(l - M R ) i R? 


e i<ut ‘ [LL-SLNLoN]! dt t 


(II.4) 


and where, 

LN = LN = LR/R 

LqN = L 0 - N (II. 5) 

— > — ) 

LL = L 0 L . 

Explicit expressions for the geometric quantities R, LN, etc. will be given below. In Eqs. 
(11.3,4) all variables subscripted "1" are evaluated at the dummy time tj; unsubscripted 
variables are evaluated at the retarded time X. 
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Given (II. 1), time can clearly be replaced by the force’s position. Thus let. 


Yi = 9-0o(M) 


Y=e- 0 o (x) 


(II. 6) 


Qt = e-e. 


With these substitutions (t! by y, , t by y, and t by 0), Eq. (II.4) becomes 


G=ie iffl(M) G(y) 

U 3 


G(y) = - 


M? LN L«N 


,icoy 


R(1-Mr) 


(II.7) 

01 . 8 ) 


T (AQ-y, ) 

+ J — n— e’^'lLL-SLN-LoNJjdy, 


(0 = co/Q 


R = ttR/U 


A0 = 0-0 O , 

Note that the retarded time relation becomes a retarded angle equation, 

A0 = y + M x R 

Substituting Eq (II.7) into Eq. (II.3) we get 


4 


PoU 3 


(H.9) 


(11.10) 


Because the right hand side of this expression depends on time only through 
0 = 0 + Qt, the time derivative can be replaced by an angle derivative (— = Q 

30 

Thus (reintroducing the pressure difference Ap) we obtain: 

fl 2 Ap iw0 o 3 


47t(uL)e‘ < ° e 


PoU 3 


— K(0-0 O , r c , r)d0 o r 0 dr 0 
30 


(11-11) 


where 
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(IU2) 


K = £ 

y 


1-M r _ 

— G . 

1-Mr 


Multiple Blades 

This gives the normal velocity induced by a single point force on one blade. If 
there are multiple (Ng) blades, let them be counted j = 1, Ng counterclockwise, so 
that 


(0 o )j = (e o )i + (H)A0b , A0g = 27t/Ng , (11.13) 

are the angular positions of all blade points at the same r 0 , x 0 . Suppose, further, that the 
loads are periodic from blade to blade, 

(Ap) j = (Ap) 1 e im(H)A0B (11.14) 

Then the net velocity induced by all Ng points is given by Eq. (II. 11) if K is 
replaced therein by 

K<e-e„| , r 0 , r) = £ e i<S " m)ae '^ l) K<e-e oj , r„, r) (11.15) 

j=l 


Geometric Relations 

For the sake of completeness we list here various purely geometric and kinematic 
relations. 

Ax = x-x 0 (t) = (yfAo) (11.16) 

where Ao = o - o Q . 

R = QR/U (11.17) 

where R 2 = (Ax) 2 + r 2 + r 2 - 2rr 0 cos y 




LN = L-N 


(11.18) 
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= [y + Ao (sin y- a cos y) - a]/R 




LqN = L 0 -N (11.20) 

= [y+ Aa - — (a*, cos y + sin y) + <Xo]/R 




LL = L 0 L 


i 1)2 

1 + —5 

Q 2 rr 0 


[(1-Kxcto) cos y+ (a-oto) sin y] 


(11.21) 




a = ra'(r), a*, = r 0 a'(r o ) (11.22) 

M r =M 0 -N (11.23) 

Q 2 rr 0 

= - M x [y + Aa + sin y]/R 


Summary 

The collective kernel function for N B blades with interblade phase index m is 
given by Eq. (11.15). The single blade kernel is obtained from Eq. (11.12) and (II. 8). 
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Appendix III 

Numerical Evaluation of the Kernel Function 


The kernel function K, defined by Eqs. (11.8), 12, 15, must be evaluated for 
specified radii (r, r 0 ) over a range of angular separations A0 max >A0> A0 m ; n . (The 
range is determined by the blade planform at radii r and r Q ). In the procedure described 
here the function K(A0) (for fixed r and r G ) is defined numerically by constructing a table 
Kj, A0j_ J = 1, N k , which spans the required domain of A0. Values of K for specific 
combinations of control and load point can then be found by interpolation. 

In the first step, the function G(y), defined by Eq. (II.8), is split into: 

G = G 0 (y) + A0(B o (y) - B 0 (— >)) - B,(y) (III.l) 


where 


M X (LN)(L 0 N) 

G 0 (Y) = 7 

R(1-M r ) 


e i<0T 


is the algebraic part of G, and: 


Y 

B„(y)= J b(Yi)/RidY, 

7m ax 

Y ^ 3 

Bi(Y)= j Yi b(Yi )/R 1 d Yi 

7m ax 

b(Y) = e'^fLL-SfLNXLoN)! . 


(III.2) 


(III.3) 

(HI. 4) 
(III.5) 


Note that "/max is an arbitrary constant and that the additive constant Bj(— °°) has been 
deleted from Eq. (III.l) (which is allowed by the fact that only changes in G are 
physically significant). 

For the reference blade the computation is begun by choosing ymax such that the 
corresponding A0 (given by Eq. (II. 9)) is just larger than the largest required value, 
A0max. The calculation is then marched backward in constant increments Ay from ymax. 
At each step the corresponding value of A0 and a provisional value of G (without the 
B 0 (— 00 ) term) are computed and saved in arrays. In addition the current value of B 0 is 
retained (so that B 0 (— °°) may eventually be found.) This procedure is terminated when 
A0 becomes less than the minimum required value. 

At this point the B 0 integration is continued (generally with a larger Ay) until the 
sending point is a prescribed axial distance downstream ( eg five tip radii.) The final 
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value of B„ is taken as B 0 (-«>) and the correction term AG B 0 (-<») is subtracted from 
each element of the G vector. The vectors A0 and G then contain, in tabular form, the 
information needed to define the influence of any point at radius r Q on any point at r for a 
given blade. 

Ifr 0 is subsonic then G is identical to the single blade kernel K (cf Eq. (11.12)). 
However, if r 0 is supersonic then A0 is a non-monotonic function of y (with extrema at 
Mr = 1) and the tables A0 and G must be interpolated onto new arrays giving K at an 
ordered sequence of A0. 

Having found K for the reference blade, the entire procedure is then repeated for 
each subsequent blade. However, the elements of the A0 vector are not the same as those 
for the reference blade. Hence for all subsequent blades the tables K, A0 must be 
interpolated onto the table for the reference blade. 

After the calculation has been completed for all blades, one has N B +1 vectors: A9 
and K for each blade. The kernel functions are then added with appropriate phasing to 
form the collective kernel, according to Eq. (11.15). 

Numerical Integration Method 

In this section the numerical method used to evaluate the integrals B 0 and Bj, 
defined in Eqs. (III.3-III.5), are described. 

Note that the integrands contain only algebraic and trigonometric functions of the 
dummy variable Yi . By integrating in equal increments. Ay, the required trigonometric 
functions can be generated recursively at each step (with some gain in computational 
efficiency). 

A simple quadrature, like Simpson’s rule, cannot be used, however, because the 

- 3 

denominator, Rj, may assume very small values over part of the integration range (when 
the control point is close to the loaded line). In order to circumvent this difficulty we use 
the rule: 


72 


I 


u(v) dy= l ay + 

R 3 (y) (R,+R 2 ) 2 -c R2 Ri 


which is exact if the numerator, u(y), is linear and R 2 (y) is quadratic, 

R 2 (y) = a + be + ce 2 , e = (y-Yi )/Ay 


(III.6) 


(III.7) 


between the limits of integration. Note that the constant and linear coefficients , a and b, 
are not needed in Eq.(III.7), but the quadratic coefficient, c, is. An expression for c will 
be worked out below. 
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Since R 2 is not quadratic, but involves a cosine of the retarded angle, the constant 
c must be determined by a curve fit. For this purpose we introduce the quadratic 
interpolation, 

cos y= (1-e 2 ) cos Yi + e 2 cos 72 + B e(l-e) (III. 8) 

e = (y-Yi )/ay> in [0,n 
B = —At[ 1 + (Ay) 2 /1 2] sin 7! . 

(The maximum error over the interval involved in this expression is (Ay ) 3 sin 7 i/72aT 3.) 
The exact expression for R is (cf Eq s . (II . 1 2, 1 3)) 

R 2 = (yfAo) 2 + D 0 - D! cosy (III-9) 

where D 0 , Di and Aa depend only on radius. 

Inserting Eq. (III. 8) into (III.9) yields the constant c required in Eq. (III. 6): 

( 

c = (A7) 2 - D] [cos 72 - cos 7 j + Ay( 1 + ^ — ) sin 7! 1 (III. 10) 

Note that c is small and so has an effect on the integral only if Rj and R2 are 
themselves small. 

Formulas III.6 and III. 10, together with a recursive evaluation of sin and cos, 
constitute the algorithm for computing B c and Bj. 
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Appendix IV 
Efficiency Calculation 


The steady state efficiency of a propeller is the ratio of the work done in forward 
motion to that done in rotation, 


UT 

UT+P l 


(IV.l) 


where P L is the power loss from induced and viscous drag ( (UT+P L )/Q is the shaft 
torque). Neglecting viscous losses, Pl is the rate of work done by the pressure forces 
acting over the blade. 


Pl = - jj pivuclA . (IV.2) 

Note that Pl is quadratic in the disturbances while UT is linear (so T) will generally 
be close to unity). In linear aerodynamic theory the errors in thrust (resulting from the 
linearization) are quadratic and therefore comparable to P L . The efficiency can be 
computed using linear theory only because ri depends on the ratio P L /UT, which is 
linear. 

In linear theory the blade is represented by a helical sheet on which there is a load 
Ap and a prescribed normal velocity. Part of the power loss, then, comes from the 
distributed loads on the sheet, 

Pls = ~jj Ap L-u r dr d0 (IV. 3) 

— ) 

where L is the normal to the helical surface. This part can be computed quite simply 

once the load distribution has been found. 

In addition to the distributed loads on the surface, however, there are concentrated 

forces acting at the leading edge, caused by the locally infinite velocities and pressures. 

These forces are quadratic in the disturbances levels and so contribute negligibly to the 

overall forces on the blade. However they act generally in the direction of motion (rather 

than at the right angles at it, like the surface forces) and so produce a contribution to the 

power which is comparable to Pls- 

— * — > 

Let U be the blade velocity and dF L an incremental leading edge force. The net 
power loss generated, then, is 

Pll = - / U dF L . (IV.4) 

LE 

Our object here is to relate this quantity to the computed load distribution, Ap, and the 
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leading edge geometry. 

Let s L (r) be the arclength along, and n^r) the unit vector normal to the leading 
edge at radius r (we take “n L to be on the helical surface and directed away from the 
blade). From a momentum balance in the locally two dimensional flow in the plane 
perpendicular to the leading edge, we find that the force on an increment of length ds^ is: 

dF L = kpPl K? ds L (IV. 5) 

where 

P L = U-(<M./ao> 2 l' 4 (IV.6) 

is the normal Prandtl-Glauert parameter. The factor K v measures the strength of the 
leading edge velocity singularity. Let p be the distance from the leading edge measured 
inboard along the local normal Iil. Then, if ul(p) is the velocity in the normal direction, 

K v = lim (ulV"|x) . (IV.7) 

|i— K) 

Note that the force dF L is in the direction of motion so that Pll is negative (i.e. the 
leading edge forces act like a power source). 

The above expression for the leading edge force is a simple restatement of Eqs. 
(9.15) and (9.16) in Ref.[12]. To be useful, though, the coefficient K v must be calculated 
from Ap. 

— > 

Let ^ be the distance measured inboard from the leading edge in the U direction 
and let A be the angle between r?L and U (i.e. A is the leading edge sweep angle). 
Referring to Fig. (1) we see that p = % cos A. 

Define a pressure coefficient singularity strength K p by 

K. = lim mT | % -^-1 . CV.8) 

5-*> p|U| 2 

From Bernoulli’s equation, then, the velocity potential behaves like 

<|>-HU|KpW (IV-9) 

Therefore, from Eqs. (IV. 8) and (IV. 9), 

K v = lim aTp-^ = ^- |U|K p a/~cosA. (IV.10) 

dp 2 y 

Finally, since = | U | cos A, we obtain the incremental power in the form, (from Eq. 
(IV. 5) and (IV. 10)) 
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(IV.ll) 


UdF L = jPP l |U| 3 k 2 ds L . 

This expression is convenient when (as here) the load Ap is computed in strips running 
— ) 

parallel to U. 

Finally we give some useful geometric results for the propeller. Recall that the 
helical surface is on x = U/Q (0 + a(r)). Let 0 l and xl be the leading edge coordinates. 
Define 


7L=r 


d0L 

dr 


S r = Q. r/U . 


Then 


(IV. 12) 


- 7 — = H + (Yl+oc) 2 /S? + y£l 1/2 


(IV. 13) 


defines distance along the leading edge, and 




|U| 

Q 


( 0 - 0 L ) 


(IV. 14) 


defines distance inboard from the leading edge (recall that a = rda/dr, as previously 
defined). 

The component of velocity normal to the leading edge is given by 

(U-n L ) 2 = ( | U | 2 + a 2 U 2 )(-^r 2 (IV. 15) 

dr 


(which can most readily be shown by first evaluating the component along the edge). 

Finally we give the two power losses in terms of the nondimensional load, 
P = S 2 Ap/pU 2 and nondimensional normal velocity W = 47tL-"u/U, which are the 
variables used in the lifting surface integral equation. 


Pls 

pQ 2 r?U 


— [ — ff PWdOrdr 
4ttS 4 l 1 


(IV. 16) 


Pll 

pQ 2 r?U 


^H PlC 


dS L 

dr 


dr 


(IV. 17) 


where 
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(IV. 18) 


c = lim [(e-ejp 2 ! 

p£ = 1 - M 2 (l + a 2 + ^S 2 )/(-^) 2 . 

The total inviscid power loss is Pl = Pll + Pls- 

In order to implement this scheme for computing efficiency the leading edge 
singularity strength C must be extracted from the numerical data. Let Pi be the value of 
P on the leading edge panel, which has a length A0i . The associated force is Pi A0) . The 
same net force is produced by the variable load 

P=yPl [ (A0O/ (0-e L ) ] 5 , (IV. 19) 

which implies a leading edge strength 

C=-i-P?A0,. (IV. 20) 

4 

Note that this is equivalent to saying that the load Pi is the load at the 1/4 chord of the 
first panel. (The same argument would assign the inboard loads to points close to the 
midchord of each panel.) The accuracy of this numerical method for computing the 
leading edge singularity strength is discussed in Appendix V. It is shown there that 
although the panel scheme converges almost everywhere, the leading edge strength has a 
residual error of about 10% in the limit of increasing panel density. It has been found that 
increasing C by a factor 1.1 significantly improves the prediction of leading edge suction, 
and therefore of induced drag. 

In summary, the efficiency is computed from Eqs. (IV. 1, 12, 13, 16, 17, 18 and 20). 
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Appendix V 
2-D Airfoil Test Case 


The constant load panel method was tested on the simple case of a 2-D thin airfoil 
in steady incompressible flow, for which an exact solution is well known. This test sheds 
some light on the selection of control point position and the identification of the leading 
edge singularity strength (which is required to calculate propeller efficiency.) 

The vertical velocity v(x) induced by a unit length vortex sheet of strength y(x) is 




(V.l) 


If, for example, we require that v(x) = 1 for x in (0,1) then the solution which obeys the 
Kutta condition y(l) = 0 is, 

7(X) = r(l-x) / Xl 5 (V.2) 

so that the total circulation and leading edge singularity strength are, 

l 

r = J y(x)dx = y (V.3) 

O 

C = V""x y(x) | x= o = 1 (V.4) 


We will examine the ability of the constant strength panel method to predict this 
solution. Using N equal length panels and a control point located a fraction e back from 
the leading edge of each panel, the corresponding discrete problem is 

-IrT ln k k ~ j+E 1 l =1 i k=l,N (V.5) 

7t rj J k-j+e-1 

r N = jv 2 if (V.6) 

W j=l 

For the leading edge singularity we examine two approximations, 

C?=y?/2>Tn (V.7) 

c? = rtf - -f- (7? + Y2 )W~ N . (V.8) 

V 8 

The first of these (which is the direct analogue of Eq. (IV. 20)) is obtained from equating 
the first panel circulation, y^Ax, to the circulation for the distribution C/Vx. Eq. (V.8) 
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is obtained by the same matching for the first two panels with an assumed distribution 
C/VW (1+bx). 

We first examine the convergence of the scheme with respect to net circulation T 
(global convergence), and the choice of control point. 

Note, first, that the solution for N = 1 is, trivially, 

T 1 =yj =7c/ln|e/(e-l)| (V.9) 

This is singular at e = 1/2 and antisymmetric in reflection e— >1— e. In fact the same thing 
is true for any value of N - so that the numerical method cannot converge to the desired 
solution for all e. 

Note, secondly, that T 1 =k/2 if e = 1/(1 ±e -2 ) = .88 or 1.16. The region around 
.88 is examined in detail in Fig. (V.l), which shows the percentage error 2 F n /tc— 1 for a 
range of N and e. The error evidently decreases to 0 as N increases for all values of e 
shown, indicating that the scheme in fact does converge to the desired solution. 

It appears that the scheme will converge (to Eq. (V.2)) for any £ > For e < 
it must converge to the image solution with Kutta condition at the leading edge. (As 
noted before the solution is always singular at e = -y .) 

Again referring to Fig. V.l, it is apparent that for each N there is a point of zero 
error, which drops from e = .88 when N = 1 to around .845 when N is large. The choice 
e = .85, which was (by a similar study) found to be optimal in the propeller problem, is 
also optimal here, producing no more than 0.1% errors in T if N > 8. 

While the scheme converges globally (whenever e > .5), it does not appear to 
converge locally at the leading edge (where the limit solution is infinite.) This is 
illustrated in Fig. V.2, which shows predictions of C from Eqs. (V.7,8) for £ = .85 and 
.82. The various curves shown do not approach the correct value (1) at large N, or even a 
common value. Furthermore, the "second order" fit, Eq. (V.8), is not noticeably more 
accurate than Eq. (V.7). 

Despite the lack of convergence, the results (for C 2 ) are within 10% of the correct 
value. This kind of error in leading edge singularity strength implies acceptable accuracy 
in the prediction of propeller efficiency. 
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Appendix VI 
Aeroelastic Analysis 

Although the aeroelastic analysis of the rotor is not strictly the topic of this report, 
the code in which the method under discussion is implemented contains options for doing 
both forced response and stability calculations. Therefore the methods used will be 
described here. The reader who wishes to supply his own aeroelastic analysis and only 
use the aerodynamic solutions can easily do so. 

We assume that the blade deformation is described by a generalized coordinate 
vector q,of length N mode , with corresponding mass, damping and stiffness matrices M,C 
and K respectively. The equations of motion of the blade are 

Mq + Cq + Kq = F(t) (VI. 1) 

where the i’th generalized force F; is the virtual work done by the aerodynamic load Ap 
against a unit amplitude deformation in the i’th mode, 

F; (t) = JJ 8;ApdA (VI. 2) 

In general we can take the force to consist of motion independent and motion 
dependent parts, 

t 

F = F 0 + j A(t-x)q(x) dx (VI.3) 

— oo 

where F 0 is the motion independent force and A is the generalized derivative of the 
force with respect to the modal coordinates ( or .equivalently, the force induced by a unit 
step in a modal displacement.) Although we cannot easily determine the detailed form of 
A from the present method, it is important that it depends only on the time lag t-T rather 
than on t and x separately. 

Forced Response 

In the forced response case we take F 0 to be simple harmonic with some prescribed 
frequency co and interblade phase index m. The solution of Eq. (VI. 1) will then have the 
same structure, so that, 

F 0 = Re(F 0 e io>l ) (VI.4) 

q = Re(qe icot ) (VI.5) 

where the complex amplitudes are related by. 
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q = [ -M co 2 + i a) C + K - Q I -1 F 0 


(VI. 6) 


The quantity Q is the matrix of frequency domain generalized forces that is computed for 
the given frequency and phase of the forcing term. Note that, from Eq.(VI.3) , the matrix 
Q is simply ico times the Laplace transform (with s = i to) of the generalized derivative 
matrix A. 

Stability 

In the stability analysis, we take Fq = 0 , and consider the temporal response to 
initial conditions q = q = 0 ,for t < 0 ; q = 0 ,q = vq for t = Of. This initial value problem 
can be solved formally by Laplace transforms (which we denote by L[] ), since the 
aerodynamic forces are in the form of a convolution. The resulting solution for 
q(s) = Lfq(t)] , at a fixed interblade phase angle, is 

q (s) = [ M s 2 + C s + K - Q(s) ] _1 M v 0 (VI. 7) 

where Q(s) = s L[A], The Laplace integral definition of Q is convergent everywhere in 
Re(s) > 0, but not in Re(s) < 0. In order for residue theory to be useful in finding the 
inverse transform of q (ie. q(t)), we must define Q(s) in Re(s) < 0 as the analytic 
continuation from the right half s plane. If this step is made, then stability can be 
determined simply by finding that root s c of the determinate, 

D(s) = det [ M s 2 + C s + K - Q(s) ] = 0 (VI.8) 

which has the largest real part. If Re(s c ) > 0 the rotor is unstable. If Re(s c ) < 0 then the 
rotor is stable for the given phase angle. In general, each interblade phase angle must be 
examined independently to find the most unstable collective mode. 

The procedure used to do the analytic continuation of Q(s) involves simple 
interpolation on the imaginary axis. The generalized forces are computed at N w user 
chosen frequencies co^ , k = 1 , N w (which typically are near one or more invacuum 
natural frequencies.) Within this frequency band Q(s) is defined by an interpolant, 

Q(s) = Ii rai GKs) ] (VI. 9) 

where the G; are analytic shape functions, and the a; are coefficients determined from the 
interpolation conditions at the knots s k = i (0* 

Li [a; G(s k) . ] = Q(s k ) (VI.10) 

Although many different shape functions could be used, the ones chosen here are, 
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G; = fl,s,s 2 , ... ,s/(s+pj_ 3 )l 


(VI.ll) 


where Pi_ 3 are real positive numbers selected to scale with the chosen frequencies. 
These functions have the property that the roots of Eq(VI.8) can be found by any 
standard complex eigensolver on an extended state space of dimension 
N s = N mo< j e max(2,N (0 -l ) consisting of the original states q , q plus ( if > 3 ) extra 
’aerodynamic states’ q ; k = 4, ... .N^. 

The eigenvalues so determined will be N s in number, and nonconjugate. Only 
those roots that fall within the circle of diameter (Dj^ - centered at i [co Nu + coj ] / 2 
are admitted. How many are found, of course, depends on how broad the original 
frequency band is. A simple strategy is to partition the range of invacuum natural 
frequencies included in the model, sweeping through the set successively with N w = 3. 
This does not increase the size of the state space for any one band, and usually leads to 
Nmode acceptable roots. If this procedure is used repetitively to find a critical ( or flutter) 
point, the frequency and interblade phase angle bands that need to be investigated can be 
narrowed very quickly as more is learned about the system. 
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Table 1 


PURDUE MODEL 


r 

A0 

275 

0.396 

3475 

0.309 

4200 

0.241 

4925 

0.180 

5650 

0.122 

6375 

0.075 

7100 

0.014 

7825 

-0.034 

8550 

-0.086 

0 

-0.163 


BLADE GEOMETRY 

b/r t \ 

.333 0. 

.333 0. 

.333 0. 

.333 0. 

.333 0. 

.333 0. 

.333 0. 

.333 0. 

.333 0. 

.333 0. 
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Table 2 


NACA 109622 GEOMETRY 


r 

AS (rad) 

b/r t 


25 

0.445 

.24 

0 . 

35 

0.315 

.24 

0 . 

45 

0.227 

.24 

0 . 

55 

0.140 

.24 

0 . 

65 

0.061 

.24 

0 . 

75 

0 

.24 

0 . 

85 

-0.052 

.24 

0 . 

95 

-0.104 

.24 

0 . 

0 

-0.131 

.24 

0 . 
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Table 3 



SR2 GEOMETRY 


r 

Ag(rad) 

b/r t 

\ 

.3 

.262 

.3 

-.036 

-4 

.215 

.3 

.091 

.5 

.150 

.298 

.150 

.6 

.087 

.296 

.173 

.7 

.024 

.292 

.145 

.8 

-.028 

.284 

.091 

.9 

-.079 

.258 

.041 

.95 

-.096 

.226 

.018 

.0 

-.122 

.102 

.009 
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Figure 2. - Helical surface construction. 


Figure 4. - Convergence of sectional thrust at 0.8 for vibrating pro- 
peller. 














Figure V.2. * Behavior of singularity strength as N increases. 


.85 .90 

E 

Figure V.1 . - Convergence of net circulation. 
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